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ABSTRACT 

The most massive stars can form via standard disk accretion - despite of the ra- 
diation pressure generated - due to the fact that the massive accretion disk yields a 
strong anisotropy in the radiation field, releasing most of the radiation pressure per- 
pendicular to the disk accretion flow. Here, we analyze the self-gravity of the forming 
circumstellar disk as the potential major driver of the angular momentum transport in 
such massive disks responsible for the high accretion rates needed for the formation of 
massive stars. For this purpose, we perform self-gravity radiation hydrodynamics sim- 
ulations of the collapse of massive pre-stellar cores. The formation and evolution of the 
resulting circumstellar disk is investigated in 1.) axially symmetric simulations using 
an a-shear-viscosity prescription and 2.) a three-dimensional simulation, in which the 
angular momentum transport is provided self-consistently by developing gravitational 
torques in the self-gravitating accretion disk. The simulation series of different strength 
of the a-viscosity shows that the accretion history of the forming star is mostly indepen- 
dent of the a-viscosity-parameter. The accretion history of the three-dimensional run 
driven by self-gravity is more time-dependent than the viscous disk evolution in axial 
symmetry. The mean accretion rate, i.e. the stellar mass growth, is nearly identical to 
the a-viscosity models. We conclude that the development of gravitational torques in 
self-gravitating disks around forming massive stars provides a self-consistent mechanism 
to efficiently transport the angular momentum to outer disk radii. Also the formation 
of the most massive stars can therefore be understood in the standard accretion disk 
scenario. 

Subject headings: Accretion disks - Instabilities - Stars: formation - Stars: massive - 
Hydrodynamics - Methods: numerical 
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1. Introduction 

The formation of circumstellar disks is a natural outcome of pre-stellar core collapse simulations 
regarding the formation of stars due to angular momentum conservation. In Kuiper et al. (2010a) 
we demonstrated the possibility, how to overcome the well-known radiation pressure barrier in the 
formation of massive stars via disk accretion. These simulations were performed in axial symmetry, 
in which the angular momentum transport in the circumstellar disk relies on a standard a-shear- 
viscosity prescription. In the three-dimensional simulation presented herein, we are now able to 
check these presumptions by studying the accretion rate driven self-consistently by gravitational 
torques, which develop from the non-axially symmetric disk structure. 

Once a nearly Keplerian disk has formed, further accretion onto the proto-star is only possible 
if angular momentum is either removed from the disk or efficiently transported to outer disk radii. 
Proto-stellar jets, outflows and disk winds allow to remove a substantial fraction of the angular mo- 
mentum from the star-disk system. Developing convective (Kley et al. 1993a; Lin et al. 1993), radi- 
ation hydrodynamical (Klahr & Bodenheimer 2003), magneto-rotational (Balbus & Hawley 1991; 
Hawlcy & Balbus 1991) and self-gravitating instabilities (Casscn et al. 1981; Anthony & Carlberg 
1988; Lin Sz Pringle 1990; Yang et al. 1991; Papaloizou & Savonije 1991; Tomley et al. 1991; Heemskerk et al. 
1992; Tomley et al. 1994; Laughlin &; Bodenheimer 1994; Miyama et al. 1994; Papaloizou &i Lin 
1995a; Laughlin & Korchagin 1996; Laughlin & Rozyczka 1996; Bate 1998; Yorke & Bodenheimer 
1999; Vorobyov & Basu 2006, 2007, 2010) in the accretion disk will transfer angular momentum to 
outer radii. Turbulent viscosity as the main mechanism for angular momentum transport in disk 
has been studied in axially symmetric simulations, including radiation transport (Rudcn Pollack 
1991; Kley & Lin 1992; Klcy et al. 1993b; Rozyczka et al. 1994), in non-axially symmetric two- 
dimensional simulations in the thin disk approximation (Vorobyov & Basu 2009) as well as in 
three-dimensional adiabatic disk simulations (Kley et al. 1993a; Lin et al. 1993). 

Self-gravity has been shown to be a major driver of angular momentum transport via develop- 
ing gravitational torques. The numerous studies on self-gravity mentioned above involve a variety 
of different approaches and methods: When starting from a collapse of a (1 to 10 Mq) pre-stellar 
core to compute the formation of the accretion disk consistently in the model, semi-analytical work 
(Lin & Pringle 1990), two-dimensional grid based hydrodynamics simulations (Laughlin & Bodenheimer 
1994; Yorke & Bodenheimer 1999), including radiation transport, and follow-up three-dimensional 
SPH disk simulations Laughlin Sz Bodenheimer (1994), as well as three-dimensional SPH simula- 
tions of the collapse (Bate 1998) have been considered. 

Several methods have been used when studying the angular momentum transport with a disk 
model as an initial condition: The first simulations performed (Cassen et al. 1981; Anthony & Carlberg 
1988; Tomley et al. 1991, 1994) were two-dimensional N-body simulations. Adams et al. (1989) and 
Laughlin & Korchagin (1996) did semi-analytical work. Papaloizou & Savonije (1991), Heemskerk et al. 
(1992), and Laughlin & Rozyczka (1996) used two-dimensional grid based hydrodynamics with a 
polytropic equation of state. The first three-dimensional grid based hydrodynamics were per- 
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formed by Yang et al. (1991), Miyama et al. (1994) performed two-dimensional SPH simulations. 
Self-gravitating disks were studied with and without the effect of an additionally small amount of 
viscosity, using a barotropic equation of state and taking into account the vertical disk structure 
(Papaloizou &: Lin 1995a). 

Several reviews outlined the importance of the ongoing research on the angular momentum 
transport in disks. The basic mechanism for angular momentum transport in disks due to gravita- 
tional torques has already been reviewed by Toomrc (1977). Tscharnutcr & Boss (199-3) overview 
theoretical core collapse models with a focus on the angular momentum transport by either turbu- 
lent viscosity (two-dimensional models) as well as gravitational torques (three-dimensional models). 
Bodenheimer (1995) present a review on observations and angular momentum transport mecha- 
nisms. Papaloizou & Lin (1995b); Lin & Papaloizou (1996) comprise angular momentum transfer 
processes in a variety of astrophysical contexts. A more recent review on angular momentum trans- 
port in accretion disks is presented by Balbus (2003) and an overview of standard accretion disk 
theory in general is given in Lodato (2008). 

The studies mentioned above focus on the collapse of a pre-stellar core similar to the solar 
nebula yielding the formation of a sun-like star. Here, we present an investigation of the angular 
momentum transport in accretion disks around massive stars. We follow the collapse of massive pre- 
stellar cores to compute consistently the formation of circumstellar accretion disks around massive 
stars. In the case of axially symmetric disks, we follow their evolution through the whole stellar 
accretion epoch. For this purpose, we make use of our recently developed self-gravity radiation 
hydrodynamics code including frequency dependent stellar luminosity feedback. The problem of 
angular momentum transport in such accretion disks around massive stars is studied both in axially 
symmetric runs by using an a-viscosity parameterization and in three-dimensional runs, in which 
gravitational torques develop from the non-axisymmetric structure and drive the accretion flow 
through the self-gravitating disk. 

Observations of this early phase of massive star formation suffer mainly from the high extinction 
of the stellar environment as well as the long distance to massive star forming regions. Recent and 
upcoming future generations of space telescopes and interferometric systems such as the Herschel 
Space Observatory, the James Webb Space Telescope (JWST) and the Atacama Large Millimeter 
Array (ALMA) will provide a deeper insight into the mechanisms of massive accretion disks. A 
list of disk candidates around massive proto-stars is already published in Cesaroni et al. (2007). 
Recently, Beuther et al. (2009) conclude from a study of twelve disk candidates that the inner 
accretion disks, whose existence are assumed due to the high coUimation degree of observed jets 
and outflows, have radii below 1000 AU and are fed by the in-falling outer envelope. So far, these 
observations fully match into the disk accretion scenario for the formation of massive stars, which 
we propose in this study. 
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2. Method 
2.1. Equations 

To follow the motion of the gas, we solve the equations of compressible self-gravity radiation 
hydrodynamics: 

dtp + V-{pu) = (1) 

dt{pu) + V -{puu + P) = f (2) 

dtE + V-iiE + P)u) = u-f-V-Ftot (3) 

V^^ = 4TTGp-2GMjr^ (4) 

dtE^ + fc^-F = (^V • i? - Q+) (5) 

F,{u,r) /F,{y,R,) = (i^./r)^ exp (-r (i/, r)) (6) 

Eqs. (1) to (3) are the conservation equations of hydrodynamics. The evolution of the gas density 
p, velocity u, pressure P, and total energy density E is computed using the open source magneto- 
hydrodynamics code Pluto3 (Mignone et al. 2007), including full tensor viscosity. The force density 
vector 

/ = -pW + Vn- AVSr- V- (F,/c)e; (7) 

includes the additionally considered physics to the equations of gas dynamics such as gravity, shear 
viscosity (in 2D only), and radiation transport. To close the system of the hydrodynamical Eqs. (1) 
to (3), we use an ideal gas equation of state 

P= (7-l)i?i„t, (8) 

which relates the gas pressure P to the internal energy E\^x_ = E — O.Bpu'^. The adiabatic index 7 
is set to 5/3. 

Eq. (4) is Poisson's equation with the gravitational potential <5 and the gravity constant G. 
denotes the central stellar mass. Our implementation of Poisson's Eq. (4) via a diffusion Ansatz 
is presented in Kuipcr ct al. (2010a). 

Eqs. (5) and (6) describe the radiation transport. In our recently developed frequency depen- 
dent radiation transport method (Kuipcr et al. 2010b) the flux of the total radiation energy density 
-Ftot is split into the flux of the frequency averaged (hereafter gray) thermal radiation energy den- 
sity F and the flux F^,{u,r) of the frequency dependent stellar irradiation. Eq. (5) denotes the 
evolution of the thermal radiation energy density E'r. The factor /c = [cyp/AaT^ -|- l) ^ depends 
only on the ratio of internal to radiation energy and contains the specific heat capacity cy and 
the radiation constant a. The source term depends on the physics included and contains any 
additional energy source such as hydrodynamical compression —PV ■ u and viscous heating. We 
solve Eq. (5) by using the so-called flux-limited diffusion approximation (hereafter ELD), in which 
the flux is set proportional to the gradient of the radiation energy density {F = —DVE^). The 
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diffusion constant D = Xc/pK^ depends on the flux limiter A, the speed of light in vacuum c, and 
the Rosseland mean opacity kr. We use the flux limiter by (Lcvcrmorc & Pomraning 1981) and 
neglect scattering. 

Eq. (6) calculates the flux of the frequency dependent stellar irradiation in a ray-tracing step. 
The flrst factor on the right hand side describes the geometrical decrease of the flux proportional 
to r~^. The second factor describes the absorption of the stellar light as a function of the optical 
depth r(i/, r) = K{v)p{r)r depending on the frequency dependent mass absorption coefficients 
For this purpose, we use tabulated dust opacities by Laor & Draine (1993), including 79 frequency 
bins, and calculate the local evaporation temperature of the dust grains by using the formula of 
Isella & Natta (2005). The flux at the inner radial boundary is given by the luminosity L*, temper- 
ature r^,, and radius i?* of the forming star. For this purpose, we use tabulated stellar evolutionary 
tracks for accreting high-mass stars, recently calculated by Hosokawa & Omukai (2009). The gas 
and dust temperature T is flnally calculated in equilibrium with the combined stellar irradiation 
and thermal radiation fleld 



Kp(r) c 



with the Planck mean opacities Kp. 



Numerical details, test cases, including a comparison of gray and frequency dependent irradi- 
ation, as well as performance studies of our recently developed hybrid radiation transport scheme 
are summarized by Kuipcr ct al. (2010b). The viscosity prescription as well as the tabulated dust 
and stellar evolution model are presented in Kuiper et al. (2010a). 

To limit the range of densities, the so-called floor value of the density is chosen to be po = 
10^21 g cm~^. This floor value occurs during the simulations only in regions where the radiation 
pressure driven outflow is depleting the density of the corresponding grid cells in the radially 
outward direction. Thus, the choice of the floor value does not influence the level of accretion onto 
the newly forming star we investigate. 

To control the spatial resolution, which is necessary to resolve the physics of self-gravity cor- 
rectly, e.g. preventing artificial fragmentation, we monitor the so-called Truelove criterion, derived 
in Truelove et al. (1997). 

In Eq. (7) the viscosity term VII is added to the conservation equations of hydrodynamics 
Eqs. (2) and (3) in case of axially symmetric, two-dimensional simulations only. The components 
of the viscous stress tensor 11 are given (in Cartesian coordinates) by 

2 
3' 

with the shear viscosity the bulk viscosity r^bi ^-nd the Kronecker symbol 6ij. We assume for the 
bulk viscosity r/b = 0. The (shear) viscosity is described via the so-called a-parameterization of 
Shakura ik Sunyaev (1973) and is computed via 

r] = pa n^ir) {H / Rf (11) 



Ilij = rj djUi + diUj - -dijdkUk + r]^5ijdkUk (10) 
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with the Keplerian angular velocity 



17k (r) 



GM{r) 



(12) 



the cylindrical radius R = r sm{9), and the pressure scale height H. M(r) denotes the mass inside 
the radius r and is given by the sum of the central stellar mass and the integral over the gas density 
in the computational domain 



In fact, the resulting rotation profiles of the forming disks in our simulations are very close to 
Keplerian {{v — vk)/vk < 2%). Of course, the actual outer edge rdisk of the accretion disk, in which 
Keplerian equilibrium is nearly satisfied, expands in time. 

If the viscosity is an effect of turbulent transport of angular momentum, it is observed that the 
strength of the stresses is proportional to the thermal pressure. This is the fundamental assumption 
of the a-Ansatz by Shakura & Sunyaev (1973). This relation holds because hotter and thicker disks 
can support higher levels of turbulence. The situation is reversed for self-gravitating disks. Here, hot 
disks are usually Toomre stable and will not produce gravito-turbulence. On the contrary, the disks 
will cool down to the marginally unstable Toomre values and establish a turbulent state where the 
level of turbulence is set by the equilibrium of energy release and radiative cooling (Gammic 2001). 
For that reason, we choose a viscosity prescription independent on the actual disk temperature 
(e.g. a fixed H/R ratio of 0.1) but only on the local mean (Keplerian) rotation profile. This way, 
we ensure that cool and thin disks can obtain the high viscosity values they deserve. Our Ansatz 
is equivalent to the so-called /3-viscosity Ansatz for self-gravitating disks by Duschl et al. (2000), 
which is also independent on temperature. A detailed derivation of the viscosity prescription is 
given in Kuipcr et al. (2010a). 



The simulations are performed on a time independent grid in spherical coordinates with a 
logarithmically stretched radial coordinate axis. The outer core radius is fixed to rmax = 0.1 pc. 
The inner core radius is fixed to rmin = 10 AU. The accurate size of this inner sink cell was 
determined in a parameter scan presented in Kuiper et al. (2010a), Sect. 5.1. The polar angle 
extents from 0° to 90° assuming midplane symmetry. In the three-dimensional run, the azimuthal 
angle covers the full domain from 0° to 360°. The grid consists of 64 x 16 x 64 grid cells, i.e. the 
highest resolution of the non-uniform grid is chosen to be 




(13) 



2.2. Numerical configuration 



Ar X rAO x rA(/>sin(6') = 1.27 x 1.04 x 1.04 AU^ 



(14) 



in the midplane {9 = 90°) around the forming massive star. The resolution decreases logarith- 
mically in the radial outward direction proportional to the radius. The radially inner and outer 
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boundary of the computational domain are semi-permeable walls, i.e. the gas is allowed to leave the 
computational domain (due to radiative forces) but cannot enter it. This outer boundary condition 
allows that the mass reservoir for stellar accretion can be controlled by the initial choice of the 
mass of the pre-stellar core. 

The Pluto code uses high-order Godunov solver methods for computing the hydrodynamics, 
i.e. it uses a shock capturing Riemann solver within a conservative finite volume scheme. The 
numerical configuration of our simulations makes use of a Strang operator splitting scheme for 
the different dimensions (Strang 1968). Our default configuration consists further of a Harten-Lax- 
Van Leer Riemann solver and a so-called "minmod" flux limiter using piecewise linear interpolation 
and a Runge-Kutta 2 time integration, also known as the predictor-corrector-method; for compar- 
ison please see van Leer (1979). Therefore the total difference scheme is accurate to second order 
in time and space. 

The internal iterations of the implicit solver for the FLD Eq. (5) is stopped at an accuracy 
of the resulting temperature distribution of AT/T < 10~^ or AT < 0.1 K. The internal iterations 
of the implicit solver for Poisson's Eq. (4) is stopped at an accuracy of the resulting gravitational 
potential of A$/$ < 10"^. 

2.3. Initial conditions 

The basic physical initial conditions are very similar to the one used by Yorke & Sonnhalter 
(2002). We start from a cold (Tq = 20 K) pre-stellar core of gas and dust. The initial dust to gas 
mass ratio is chosen to be Mdust/-^gas = 1%- The model describes a so-called quiescent collapse 
scenario without initial turbulent motion (ur = ug = 0). The core is initially in slow solid-body 
rotation (|ti^|/i? = Q.q = 5 * 10"^'^ Hz). In the axially symmetric runs, the total mass A^core i^ 
the computational domain is 60 or 120 Mq (17.2 or 48.6 Mj, respectively) and the value for the 
physical a- viscosity varies from a = up to a = 1. In the three-dimensional run, the total mass 
Mcore in the computational domain is chosen to be 120 Mq (48.6 Mj). An overview of the runs 
performed is given in Table 1. For a more comprehensive parameter scan of the initial conditions 
of the axially symmetric core collapse simulations, please see Kuipcr ct al. (2010a). 
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Table 1. Overview of simulations presented. 



Label 


Dimension 


-^'^-'core L WJ 


(X 


'^end v'^j ^ J 


-^•^-^^V ena7 L 0J 


2D-60Msol-alphal.O 


2D 


60 


1.0 


449.0 


28.7 


2D-60Msol-alpha0.3 


2D 


60 


0.3 


939** 


28.2 


2D-60Msol-alpha0.1 


2D 


60 


0.1 


240.7 


21.7 


2D-60Msol-alpha0.05 


2D 


60 


0.05 


126.7 


16.8 


2D-60Msol-alpha0.03 


2D 


60 


0.03 


48.0 


12.6 


2D-60Msol-alpha0.02 


2D 


60 


0.02 


48.3 


12.1 


2D-60Msol-alpha0.01 


2D 


60 


0.01 


53.4 


11.4 


2D-60Msol-alphaO 


2D 


60 





24.0 


8.8 


2D-120Msol-alphal.O 


2D 


120 


1.0 


14.5 


23.2 


2D-120Msol-alpha0.3 


2D 


120 


0.3 


489.0** 


56.5 


2D-120Msol-alpha0.1 


2D 


120 


0.1 


38.7 


25.6 


2D-120Msol-alpha0.03 


2D 


120 


0.03 


34.6 


24.7 


2D-120Msol-alpha0.01 


2D 


120 


0.01 


11.3 


18.7 


2D-120Msol-alpha0.003 


2D 


120 


0.003 


10.8 


18.2 


2D-120Msol-alpha0.001 


2D 


120 


0.001 


9.3 


18.1 


2D-120Msol-alphaO 


2D 


120 





8.2 


18.0 


3D-120Msol 


3D 


120 





11.8 


24.5 



Note. — The runs differ in their dimension, their initial pre-stellar core mass Mcore; and their 
value of the a- viscosity-parameter. The period tend of evolution simulated and the corresponding 
stellar mass at tend are given. A ** denotes that the computation has been stopped at the 
point in time when no mass is left in the computational domain. 
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3. Results 

3.1. Viscous disks (2D) 

Developing convective (Kley et al. 1993a; Lin et al. 1993), radiation hydrodynamical (Klahr & Bodenheimer 
2003), magneto-rotational (Balbus & Hawley 1991; Hawley & Balbus 1991; Flock et al. 2010; Dzyurkevich et al. 
2010) and self-gravitating instabilities (Yang et al. 1991; Laughlin &: Bodenheimer 1994; Bodenheimer 
1995; Vorobyov & Basu 2006, 2007, 2010) in the accretion disk will transfer angular momentum 
to outer radii. To mimic the effect of this angular momentum transport, we made use of the a- 
viscosity model by Shakura & Sunyaev (1973). The differential rotation in the circumstellar disk 
yields a shear viscosity, transporting angular momentum to outer disk radii, while enhancing the 
stellar accretion rate. 

To study the influence of the value of the a-shear-parameter on the resulting stellar accretion 
rate, we first perform several simulations of an axially symmetric 60 Mq pre-stellar core collapse 
with varying normalization values for the physical a-viscosity. The resulting stellar mass and the 
accretion rate as a function of time is shown in Fig. 1 . 
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Fig. 1. — Stellar mass M^, (upper panel) and accretion rate M* (lower panel) as a function of time 
t for different values of the a-shear parameter during the collapse of an axially symmetric 60 Mq 
pre-stellar core. 

Remark: The a- values mentioned here are associated with a constant aspect ratio H/R = 0.1 of 
the circumstellar disk as described in Sect. 2.1. Essentially, the H/R ratio computed from the 
simulation data depends on the radius and increases slightly from 0.05 at the inner disk rim up to 
0.15 at the outer disk radius at r 200 AU. 



3.2. Self-gravitating disks (3D) 



To achieve a more detailed picture of massive accretion disks and to cross-check the assumptions 
made by applying an axially symmetric viscous disk model, we expand the two-dimensional setup 
into three dimensions now. In the three-dimensional run no shear viscosity term is applied (a = 0). 
Instead, the forming massive accretion disk will self-consistently evolve non-axisymmetric structures 
yielding gravitational torques, which efficiently transport angular momentum to outer disk radii, 
allowing further accretion onto the new born star. We perform the 3D simulation with an initial 
core mass of 120 Mq. A comparison of the resulting stellar mass growth and the accretion history 
for the 3D case as well as several 2D cases with and without viscosity is presented in Fig. 2. Due 
to the high computational costs of the radiative 3D hydrodynamics simulations from the pre-stellar 
core collapse scale down to the 1 AU scale of the accretion disk formation and evolution, we are 
currently not able to report on the long-term evolution of the pre-stellar core collapse. The 11.8 kyr 
of collapse evolution simulated correspond to 0.25 tg of the initial pre-stellar core. For comparison, 
in our previous 2D collapse simulations in Kuipcr ct al. (2010a), also the case a = 0.3 herein, we 
followed the collapse over 10 until no mass is left in the computational domain; it is either 
accreted onto the central star or expelled over the outer boundary at 0.1 pc by radiative forces. 
Here, we will focus on the formation of the disk and the evolution of the main disk accretion phase, 
which is well resolved in the 3D simulation as well. The disk accretion epoch in the 120 Mq run 
starts roughly at 7.8 kyr and we follow the evolution up to 11.8 kyr. These 4 kyr correspond to 
566 inner disk orbits (at 10 AU) or 28 outer disk orbits (on average roughly at 75 AU) respectively, 
based on a mean stellar mass of M*(7.8 — 11.8 kyr) ~ 20 Mq. 

Compared to the axially symmetric simulations of the viscous accretion disk, the angular 
momentum transport by gravitational torques due to the 3D non-axially symmetric disk morphology 
results in a much more time dependent accretion history (lower panel of Fig. 2), rather successive 
accretion events than a smooth accretion flow. On the other hand, the integral over the accretion 
history Jq M(t) di, the actual stellar mass M*(t), remains qualitatively quite similar so far and 
finally yields a slightly higher total stellar mass (upper panel of Fig. 2). 
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Fig. 2. — Stellar mass (upper panel) and accretion rate M^, (lower panel) as a function of time t 
or number of inner orbits during the collapse of a 120 Mq pre-stellar core. The angular momentum 
transport in the two-dimensional axially symmetric run (lines with symbols) is driven by shear 
viscosity, calculated via an a-parameterization prescription. The angular momentum transport in 
the three-dimensional run (solid line) is driven by gravitational torques in the self-gravitating disk, 
see Fig. 6 for a visualization of the non-axisymmetric disk structure. To visualize in particular 
the details during the disk accretion epoch, the lower panel displays the section from 7.8 kyr up 
to 12.1 kyr; as in the 60 Mq case, the accretion rate during the free fall epoch {t < 7.8 kyr) is 
constant. 
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4. Discussion 

4.1. From the free fall phase to the disk accretion epoch 

A characteristic feature of the accretion histories, see lower panel of Fig. 1, is a sharp drop 
down at roughly t = 8 kyr. This point in time marks the transition between the initial 'free fall' 
phase and the disk formation epoch afterwards. During the free fall phase, the gas dynamics are 
dominated by the gravitational drag towards the core center. During the course of the collapse, 
high angular momentum gas from initially larger radii reach the inner core. At a specific point in 
time, the resulting centrifugal force (with slight support by thermal pressure) completely balances 
the gravity at the inner computational boundary (in the first grid cell adjacent to the central sink 
cell). Hence the accretion rate drops down. 

From now on, the follow-up mass spiraling inwards builds up a circumstellar disk. In other 
words: The so-called centrifugal radius is increasing in time and therefore a circumstellar disk 
grows during the collapse from the inside outwards. Further accretion through these disks has to 
be attended by an efficient angular momentum transport. 

4.2. Viscous disks (2D) 

Our simulation series in axial symmetry with varying values of the a-shear-viscosity can be 
separated into two distinguishable regimes. The runs with low viscosity (a < 0.05) yield the 
formation of rings instead of the formation of an extended accretion disk. As stated by Yorke et al. 
(1995), these rings would break up in corresponding three-dimensional studies, resulting in angular 
momentum transport by tidal torques. Runs with higher a-viscosity yield the formation of long- 
living accretion disks. Compared to analytical estimates on the a- values of massive accretion disks 
considering typical high mass accretion rates in standard accretion disk models (including the 
thin disk approximation), presented in Vaidya et al. (2009), our dynamical study allows also the 
formation of massive accretion disks with slightly larger a- values. 

During the long-term evolution of these viscous circumstellar disks, the accretion rate onto 
the central star decreases continuously. The long-term runs differ only about 16% in the resulting 
stellar mass after 10^ yrs of evolution, although the a-values of these runs vary up to a factor of 
20. After 10^ yrs the star has grown up to 19.0, 19.1, 17.8, and 16.0 Mq for a = 1.0, 0.3, 0.1, and 
0.05 respectively. If we infer from these simulations that the accretion rate is rather independent 
on a, the well-known relation for the accretion rate in analytical steady state models 

= SttuJ: (15) 

implies that the surface density T, is inversely proportional to the dynamical viscosity u, as it is 
expected from self-similar solutions for cases, in which the viscosity has a power-law dependence 
on the radius (see e.g. Hartmann 1998). To backup this conclusion. Fig. 3 shows the mass in the 
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midplane (up to the first pressure scale height H = O.li? of the disk or /S.6 = 5.625° respectively) 
as a function of the radius for varying values of the a-parameter for two different points in time. 
Due to the fact that a is only the normalization of the strength of the shear-viscosity, the slope 
of the in-falling envelope at radii larger than the disk radius, where the shear force is negligible, is 
independent of the a-parameter. Going to smaller radii, we observe a characteristic steep increase in 
mass. This accumulation of the in-falling mass in the midplane is simply a results of the conservation 
of angular momentum. Each volume element at the initial radius r and the polar angle 9 can be 
related to its so-called centrifugal radius -Rcent? where gravity is completely compensated by the 
centrifugal force: 

^^'^^^ - G Mir) ^^^^ 

with the initial angular velocity f^o the enclosed mass M{r) inside the radius r. This equilibrium 
of gravity and centrifugal force occurs at the inner computational boundary (at 10 AU) at roughly 
7.8 kyr, afterwards the region of equilibrium increases radially. If we define the location of the 
steep increase in mass as the outer disk radius -Rdisk) we can infer a rough estimate on the growth 
rate of the circumstellar disk: In the case of a = 0.3, the outer disk radius i?disk follows up to the 
first free fall time (tg = 67.6 kyr) of the pre-stellar core roughly the linear relation 

29 AU 

^disk ~ ~r~i ~ ionsct) + -Ronsot (17) 

1 kyr 

with the beginning of the disk formation at ionset = 7.8 kyr and -Ronsct = 10 AU. It has to be kept in 
mind that the outer radius i?disk of this accretion disk is defined via the steep increase in mass and 
the outer parts of this region r < -Rdisk are not necessarily in Keplerian rotation, see Fig. 4. At the 
outer disk radius, the in- falling new disk material is not in gravitational-centrifugal equilibrium. 
The material in motion falls down to a radius smaller than its centrifugal radius and therefore this 
region of disk growth is characterized by super-Keplerian motion. For an observational example of 
super-Keplerian motion around young forming high-mass stars, see Bcuthcr & Walsh (2008). At 
larger radii, the envelope - still feeding the accretion disk - shows a strong sub-Keplerian profile. 
This outcome supports the speculation by Beuther et al. (2009), who state that no massive accre- 
tion disk in Keplerian motion could be observed so far due to the current observational resolution 
limit of roughly 1000 AU. But embedded in the currently observed large-scale flattened structures 
in non-Keplerian motion, they expect to find Keplerian massive accretion disks with radii smaller 
than 1000 AU. The size of the disks simulated here expands in time. 

In the sheared disk region (r < -Rdisk) higher values of a indeed result in lower mass disks. As 
discussed previously, for low viscosity runs (a < 0.05) the evolution leads to the formation of rings 
instead of a stable accretion disk. 

The accretion rates as a function of the radius are displayed in Fig. 5. In the low-a-simulations 
with a < 0.05, in which rings have formed rather than accretion disks, the accretion rate varies 
strongly with the radius. If the a- viscosity is high enough to maintain the accretion flow, the accre- 
tion rates are indeed approximately constant throughout the whole disk region and also independent 
of the strength of the a-viscosity. 
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Fig. 3. — Gas mass M in the domain's midplane (from midplane up to the first pressure scale height 
H = O.li?) as a function of the radius r for varying values of the a-parameter for two different 
points in time. 
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Fig. 4. — Deviation of the azimuthal velocity from Keplerian motion in the domain's midplane 
for the case of Mcore = 60 Mq and a = 0.3 for different points in time. The Keplerian velocity 
is hereby calculated based on the stellar mass plus the enclosed mass M(r) in the computational 
domain. The inner accretion disk is characterized by Keplerian motion. The outer edge of the 
growing accretion disk is characterized by a super-Keplerian peak due to the newly in-falling mate- 
rial, which is not in gravitational-centrifugal equilibrium yet. The outer envelope is characterized 
by a sub-Keplerian profile towards the initial solid body rotation (t = 0). 
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Fig. 5. — Accretion rate M integrated over the polar angle as a function of the radius r for 
varying values of the a-parameter for two different points in time. 
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In terms of the pre-stellar core collapse simulations presented here this means that the impact of 
the choice of the a- value on the accretion rate is almost negligible (unless the a-value is chosen high 
enough), but strongly affects the morphology of the disk. In the following section, we compare these 
axially symmetric a-disk models with a self-consistent scenario, in which the angular momentum 
transport is calculated from the disk evolution itself. 
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4.3. Self-gravitating disks (3D) 

Hence we started to expand our simulations into three dimensions. These simulations al- 
low to compute the angular momentum transport by developing gravitational torques in the self- 
gravitating disk consistently with the formation and evolution of the accretion disk. A visualization 
of such a non-axisymmetric density structure yielding a gravitational torque is presented for the 
case of a 120 Mq collapse in Fig. 6. The figure shows a face-on view of a slice through the inner 
core region at t = 10 kyr. 

The similarity of the mean accretion rate with the high-ce-disks justifies the assumption made 
in the viscous disk simulations, in which the computation of the azimuthal structure of the accretion 
disks is replaced by an a priori a-parameterization. The a-prescription comprises to a reasonable 
accuracy the overall angular momentum transport by evolving instabilities in the three-dimensional 
disk. 

The origin of these instabilities lies in the self-gravity of the circumstellar disk as supported 
by the fact that the Toomre Q value 

■3 = 2^ <^«' 

approaches unity from above (Q — > 1 and Q > 1) at the unstable radius during the formation of 
the spiral pattern due to an enhancement of the surface density. Fragmentation of the disk during 
the further evolution of this instability is likely to be suppressed by the heat as well as the fast 
rotation of the inner disk. The Toomre Q value in the unstable regime stays above unity so far 
(1< Q < 2). 

Modes with lower azimuthal wave number (especially m = 1) are potentially susceptible 
for a growth by the so-called SLING instability ("Stimulation by the Long-range Interaction of 
Newtonian Gravity"), which was first studied numerically in Adams ct al. (1989) and analytically 
in Shu ct al. (1990). Nonetheless, amplifications by the SLING instability require a disk to stellar 
mass ratio of Mdisk/^^star > 1/3- The SLING instability seems to be the major driver of ampli- 
fication up to the point of binary formation if Mjisk/Mstar 1) which could be the case during 
the early onset of a pre-stellar core collapse before the central star dominates (Adams & Benz 
1992). Note, that m = 1 modes cause the center of mass to be unequal to the geometrical cen- 
ter of the disk, which is not taken into account in our simulations. Anyway, in our simulations 
the instability occurs when the star has already grown up to 20 M© so that amplification by the 
SLING instability seems to be unlikely. Furthermore, the amplification by the SLING instability 
only grows slowly on timescales of a few outer rotation periods of the disk (with a minimum of 
one outer rotation period for Mdisk/A^star = !)• Such a large mass ratio is observationally detected 
when comparing the large scale (non-Keplerian) fiattened structure with the proto-stellar mass 
(e.g. Ccsaroni ct al. 1997; Zhang ct al. 1998; Bcuthcr ct al. 2005). Contrary to this slow growth, 
the so-called SWING amplification of modes with larger azimuthal wave numbers m > 1, proposed 
first by Goldreich &: Lynden-Bell (1965) and Toomre (1981), potentially will grow much faster. Es- 
pecially once axisymmetric disturbances due to self-gravity become prominent, modes with m > 1 
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Fig. 6. — Face-on view of the high density gas in a shce through the x-y midplane. The length of 
the chpping domains goes from 50 AU x 50 AU up to 400 AU x 400 AU. The data is taken from a 
snapshot {t = 10 kyr) of the 120 Mq pre-stehar core cohapse. A linear interpolation was applied to 
the discretized data cube. The density of the central sink cell with a radius of 10 AU is calculated 
by an extrapolation of the gas density at the boundary assuming a zero gradient in the density, so 
this region should be taken with care. 
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dominate (cp. Papaloizoii & Savonije 1991). 

The Fourier analysis of the non-axisymmetric disk formed in the 120 Mq collapse is displayed 
in Fig. 7. The figure shows the amplitudes of the modes for different radii and the mean amplitude 
of the sum over all radii (upper panel) as well as the evolution of the mean amplitude with time 
(lower panel). The upper panel documents that the m = 2 mode has a higher amplitude than the 
m = 1 mode at the inner rim of the disk {{A[m = 2) — A{m = \))/A{m = 1) = 26%). Elsewhere 
m = 1 is slightly higher. The evolution of the mean amplitude (normalized sum over all radii) with 
time, shown in the lower panel, demonstrates that during the onset of the gravitational instability 
the m = 2 mode grows faster than the others. The mean growth rate is roughly independent on 
the azimuthal wave number m. 

Further well-founded conclusions on the evolution (the growth of unstable modes with different 
azimuthal wave numbers) and fate (potential fragmentation) of the spiral patterns require to study 
the disk's evolution in the three-dimensional run for a much longer period (potentially up to a 
factor of 40 for a complete coverage of the stellar accretion phase) and to analyze a variety of 
initial conditions. Moreover, a study of the potential fragmentation of the three-dimensional disk 
structure is currently inhibited in these simulations due to the fact that after the onset of spiral 
arm formation the Truelove criterion is violated. The Jeans number J = Ax/Aj in the high-density 
regions approaches unity, i.e. the Jeans length Aj = ^Jirc^/Gp is only resolved by a single grid 
cell with the resolution Ax. Truelove et al. (1997) recommended a Jeans number of J < 0.25, 
i.e. to resolve the Jeans length by at least four grid cells. Hence, the best and straight forward 
improvement would be to go to an even higher resolution of the forming circumstellar disk. An 
alternative approach, used e.g. in Krumholz et al. (2007, 2009), would be to use so-called sink 
particles when the Truelove criterion is violated. In the sink particle method, the Truelove criterion 
implies a maximum density, which is considered in the gas dynamics (regions with higher density 
are represented by the particles). The method of sink particles should be justified, if it is guaranteed 
that the gas accumulated in the sink would be gravitationally bound anyway. On the other hand, 
this study focuses on the resulting accretion history of the forming massive star, which initial phase 
shown here completely relies on the well resolved large scale structure of the non-axisymmetric 
spiral arms. Fragmentation on larger scales potentially limits the mass reservoir of the forming star 
in the center (Krumholz et al. 2009; Peters et al. 2010b, c), although Peters et al. (2010a) recently 
have shown that a large scale magnetic field will offer locally support against gravitational collapse 
reducing secondary fragmentation. 

The results of our disk accretion simulations match the previous models of viscous and self- 
gravitating disks by Vorobyov & Basu (2009), proposing viscous-dominated disks for a > 10~^ and 
very high disk accretion rates for a > 1Q~^. The episodic accretion events in our self-gravitating 
three-dimensional simulations and the non-episodic accretion of our two-dimensional viscous disks 
are in good agreement with the results by Vorobyov &: Basu (2010), who state that "the (additional) 
effect of a generic a-type viscosity acts to reduce burst frequency and accretion variability, and is 
likely to not be viable for values of a significantly greater than 0.01". For a comparison of 
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Fig. 7. — Fourier analysis of the non-axisymmetric accretion disk formed during the collapse of a 
120 Mq pre-stellar core. The upper panel displays the strength of the modes m at a specific point 
in time t = 11 kyr. The different lines clarify the dependency on the radius r and denote the modes 
at the inner disk rim (dotted line with stars), at the outer core region (dashed dotted line with 
squares) as well as the mean amplitude of the total sum over all radii (solid line with diamonds). 
The lower panel shows the evolution of the mean amplitude with time for the four highest modes 
m = 1, 2, 3, and 4. The amplitudes in both panels are normalized to the m = mode. 
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these outcomes with previous semi-analytic predictions including disk fragmentation, please see 
also Krattcr & Matzner (2006); Krattcr ct al. (2008). 

A direct comparison of the resulting accretion histories to observations suffer from the fact that 
currently the stellar accretion rate cannot be observed in the case of high-mass star formation due to 
the high extinction of the proto-stellar environment. Observations with reference to low-mass star 
formation clearly favor the episodic accretion events shown in our three-dimensional simulation, 
see e.g. Enoch et al. (2009); van Boekel et al. (2010) for most recent examples. Previous three- 
dimensional simulations by Krumholz et al. (2007, 2009) fully support this outcome, too. 
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5. On the stability of radiation pressure driven outflows 

Although this pubhcation highly focuses on the angular momentum transport and accretion 
flow through the circumstellar disk, the three-dimensional simulation reveal more interesting phys- 
ical processes. 

5.1. The issue 

Background of this section is that our simulation series in axial symmetry presented in Kuiper et al. 
(2010a) have shown the launch of a stable radiation pressure driven outflow in the bipolar direction. 
As a result, the radiation pressure removes a substantial fraction (more than 40% during the stellar 
disk accretion phase) of the initial mass of the pre-stellar core from the star-disk system, reducing 
the star formation efficiency of the massive star forming core. 

In the simulations by Krumliolz ct al. (2009), these outflows are affected by the so-called '3D 
radiative Rayleigh- Taylor instability'. In our three-dimensional simulation presented here, the 
strong radiation pressure in the bipolar direction stops the in-fall and reverts the mass flow into 
a radiation pressure driven outflow with velocities of a few 100 km s~^. Fig. 8 visualizes a slice 
through the three-dimensional density structure of such an outflow, launched during the collapse 
of a 120 Mq pre-stellar core. In contrast to Krumholz et al. (2009), the outflow in our three- 
dimensional simulation remains stable despite the strong non-axially symmetric features developing 
in the outflow region so far. 

5.2. Possible explanations 

The reason for the disagreement on the stability of the radiation pressure driven outflow can 
be manifold. Although we cannot guarantee the completeness of the following listing, we try our 
best to objectively overview the most relevant explanations here: 

1. We use an initially steeper density profile (p oc r~^) than Krumholz et al. (2009) [p oc r~^'^). 
I.e., at the point in time, when the radiation pressure driven outflow is launched, there is 
more mass in the bipolar direction, which is swept up at the border of the developing cavity, 
in the configuration by Knmiholz ct al. (2009). 

2. The instability potentially requires non-axial symmetric modes (3D). Therefore our 2D runs, 
including the whole stellar accretion phase, show a stable and long-living outflow. In our 3D 
run, the outflow is stable during the launching and the first growth phase, but the instability 
will occur at a later epoch, e.g. the instability in Krumholz et al. (2009) occurs when the 
outflow extents up to slightly more than r ~ 1000 AU and the outflows in our runs have just 
propagated up to r < 1000 AU so far. 




Fig. 8. — Edge-on views of the low density gas in a slice through the x-z plane. The data is taken 
from a snapshot (t = 10 kyr) of the 120 Mq pre-stellar core collapse. A third order interpolation 
was applied to the discretized data cube. 
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3. Our simulations account for frequency dependent stellar irradiation instead of the frequency 
averaged (gray) approximation used in Krumholz ct al. (2009). In case of gray radiation 
transport, the stellar irradiation is absorbed in a thin layer (the 'radiation bubble wall' in 
Krumholz et al. (2009)) independent of its natural frequency. In case of frequency dependent 
irradiation, the IR part of the stellar spectrum is absorbed behind the UV absorption layer. 
This potentially results in a pre-acceleration of the layers on top of the actual cavity. The 
sharp wall as seen in Krumholz et al. (2009) changes to a broader transition region, if the 
frequency dependence is considered. 

4. Additionally to the frequency dependence, our simulations use a more accurate treatment 
(ray-tracing) of the 'first absorption' of the stellar irradiation in the cavity wall, and the 
instability is an artifact of the flux-limited diffusion approximation. The diffusion approxi- 
mation, the assumption that photons behave like a diffuse gas, is only valid in an optically 
thick medium. The flux limiter guarantees also for an optically thin medium that the radia- 
tion velocity is always less than or maximal equal to the speed of light, but does not represent 
the correct motion of photons: In the optically thin cavity, photons will propagate on straight 
lines impinging on the cavity wall, instead of a diffuse expansion such as a gaseous bubble. 
Exactly in the cavity wall, at r ~ 1, the flux- limited diffusion approximation breaks down 
and yields the wrong radiative flux. 

5. The disagreement is due to the different numerical resolution of the diverse grid structures: 
Under the assumption that the refinement criteria of the AMR grid guarantees the usage of 
the highest resolution level at the cavity walls {Ax = 10 AU in the case of Krumholz ct al. 
(2009)), our fixed and radially stretched grid in spherical coordinates has the same resolution 
of Ax ~ 10 AU roughly at a radius of r = 100 AU. In the inner core region (r < 100 AU) 
the resolution of our fixed grid increases gradually up to Ax r-^ 1 AU. In the outer core region 
(r > 100 AU) the resolution of our fixed grid decreases gradually down to Ax ^ 2000 AU. 
Significant statements on the resolution needed require more resolution tests with both grid 
structures. Furthermore, for a comparison of the different resolution, the following item has 
to be considered, too. 

6. We use a grid in spherical coordinates, while Krumholz et al. (2009) use Cartesian. Both 
dominant forces in the bipolar region (gravity and radiation pressure) are in first order aligned 
with the radial coordinate axis, i.e. they are better represented in a spherical coordinate 
system. This difference in the coordinate system gets unimportant by the time the resolution 
needed is known and used. But at least this points out that a higher resolution is needed when 
using Cartesian coordinates to represent the radially acting forces sufficiently (e.g. analogously 
much higher resolution is needed in Cartesian coordinates to represent a circular orbit than 
in spherical coordinates, where the orbital motion is aligned with the azimuthal axis). 



5.3. Conclusion 



The simulations by Krumholz ct al. (2009) and ourselves (Kuipcr ct al. 2010a, and the sim- 
ulations herein) differ in a large variety of numerical methods, initial conditions, and radiative 
physics included, which could be responsible for the disagreement on the stability of the outflow 
region. We suppose that the underlying physical difference is the fact that the expansion of the 
radiation bubble is stopped in the simulations by Krumholz ct al. (2009). I.e., the gravitational 
force is roughly as high as the radiation pressure in this situation (during the launch of the bubble 
the radiation pressure has to be higher). Then the marginal stable bubble wall is subject to a ra- 
diative Ray leigh- Taylor instability. But, in our opinion, this kind of instability would work in two 
dimensions as well and will not require non-axially symmetric modes to occur. In our simulations, 
the central star still accretes mass after the launch of the outflow, yielding a continuous growth 
of its luminosity. The radiation pressure remains stronger than gravity and the expansion of the 
optically thin cavity region is not stopped. In this way, the radiation pressure expels more than 
half of the initial core mass from the star-disk system and therefore efficiently reduces the star 
formation efficiency. 

To conclude: The details of the underlying physics of the '3D radiative Rayleigh- Taylor in- 
stability' are quite unknown so far and deserve further investigation. A well-founded answer to 
the stability will have an impact on the more general question of the importance of feedback pro- 
cesses on the star formation efficiency. Anyway, the final stellar masses of our axially symmetric 
simulations, including a stable long-living bipolar outfiow, evidently support the statement that 
instabilities in the outflow regions are not required to form the most massive stars. 
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6. Summary 

We performed axially symmetric and three-dimensional high-resolution radiation hydrody- 
namics simulations of monolithic pre-stellar core collapses towards the formation and evolution 
of massive accretion disks including frequency dependent radiative feedback. The vicinity of the 
forming star could be resolved down to 1.27 AU. We directly compared the history of accretion 
rates driven either by gravitational torques due to the self-gravity of the massive circumstellar disk 
or by an ct-shear-viscosity prescription. The influence of the choice of the a-value on the resulting 
accretion history was determined in a simulation series. 

We found that the accretion rates and therefore also the final stellar masses depend only 
marginally on the strength of the a- viscosity as long as a sufficient transport of angular momentum 
is guaranteed to avoid the formation of ring instabilities. Instead, the choice of the strength of 
the a-shear-viscosity determines the morphology of the forming circumstellar disk in such a way 
that higher a- values lead to lower surface densities. Contrary to the continuous accretion during 
the viscous disk evolution, the three-dimensional simulation results in time dependent episodically 
accretion events, as already observed in the case of low-mass star formation. Nevertheless, the 
mean accretion rate is very similar to the a- viscosity approach of the axially symmetric simulations 
(even a slightly higher mean accretion rate can be obtained by the developing gravitational torques) 
and therefore fully justifies the assumptions made in the two-dimensional studies. 

Our simulations show the launch of radiation pressure driven outflows, which as long as we 
can proceed with the three-dimensional run maintain their stability in both axial symmetric and 
non-axial symmetric stellar environments. 

Finally, the angular momentum transport via gravitational torques is shown to be sufficiently 
strong to allow for the high accretion rates needed during the formation of the most massive stars. 
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